home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / LUDCMP.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  65 lines
  1. PROCEDURE ludcmp(VAR a: glnpbynp; n,np: integer;
  2.        VAR indx: glindx; VAR d: real);
  3. (* Programs using LUDCMP must define the types
  4. TYPE
  5.    glnpbynp = ARRAY [1..np,1..np] OF real;
  6.    glnarray = ARRAY [1..n] OF real;
  7.    glindx = ARRAY [1..n] OF integer;
  8. in the main routine. *)
  9. CONST
  10.    tiny=1.0e-20;
  11. VAR
  12.    k,j,imax,i: integer;
  13.    sum,dum,big: real;
  14.    vv: glnarray;
  15. BEGIN
  16.    d := 1.0;
  17.    FOR i := 1 TO n DO BEGIN
  18.       big := 0.0;
  19.       FOR j := 1 TO n DO IF (abs(a[i,j]) > big) THEN big := abs(a[i,j]);
  20.       IF (big = 0.0) THEN BEGIN
  21.          writeln('pause in LUDCMP - singular matrix'); readln
  22.       END;
  23.       vv[i] := 1.0/big
  24.    END;
  25.    FOR j := 1 TO n DO BEGIN
  26.       FOR i := 1 TO j-1 DO BEGIN
  27.          sum := a[i,j];
  28.          FOR k := 1 TO i-1 DO BEGIN
  29.             sum := sum-a[i,k]*a[k,j]
  30.          END;
  31.          a[i,j] := sum
  32.       END;
  33.       big := 0.0;
  34.       FOR i := j TO n DO BEGIN
  35.          sum := a[i,j];
  36.          FOR k := 1 TO j-1 DO BEGIN
  37.             sum := sum-a[i,k]*a[k,j]
  38.          END;
  39.          a[i,j] := sum;
  40.          dum := vv[i]*abs(sum);
  41.          IF (dum > big) THEN BEGIN
  42.             big := dum;
  43.             imax := i
  44.          END
  45.       END;
  46.       IF (j <> imax) THEN BEGIN
  47.          FOR k := 1 TO n DO BEGIN
  48.             dum := a[imax,k];
  49.             a[imax,k] := a[j,k];
  50.             a[j,k] := dum
  51.          END;
  52.          d := -d;
  53.          vv[imax] := vv[j]
  54.       END;
  55.       indx[j] := imax;
  56.       IF (a[j,j] = 0.0) THEN a[j,j] := tiny;
  57.       IF (j <> n) THEN BEGIN
  58.          dum := 1.0/a[j,j];
  59.          FOR i := j+1 TO n DO BEGIN
  60.             a[i,j] := a[i,j]*dum
  61.          END
  62.       END
  63.    END;
  64. END;
  65.